function [y,K,B,idxComp,idxBn,Kbc,retBnval]=solvefem2dvec(nodesxyz,elemnode,boundrnode)
%-------------------------------------------------------------
% solvefem2dvec version 1.0, based on poissonfem2d version 1.0
% this function solves 2D electrostatic 
% it follows the algorithm provided in (Jianming Jin) 
% "The finite element method in electromagnetics"
%--------------------------------------------------------------
% Copyright 2010 Oka Kurniawan (kurniawano@ihpc.a-star.edu.sg)
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.
%-------------------------------------------------------------
% output matrices: y, K, B
% this function solves [K]{y}={B}
% input matrices: nodesxyz,elemnode,boundrnode
% => nodesxyz: contains the (x,y) positions of the nodes, 
%                      size=(number of nodes) X 2
% ========> format: column 1 = x position
% ========>             column 2 = y position
% => elemnode: contains the global index of the nodes in elements, 
%                       and element parameters
%                       size=(number of elements) X 7
% ========> format: column 1 = node 1
% ========> format: column 2 = node 2
% ========> format: column 3 = node 3
% ========> format: column 4 = alpha_x
% ========> format: column 5 = alpha_y
% ========> format: column 6 = beta
% ========> format: column 7 = f
% where alpha_x, alpha_y,beta, and f is defined in (Jianming Jin)
% => boundrnode: contains nodes for the Dirichlet Boundary conditions,
%                          and its values, 
%                          size = (number of nodes in dirichlet b.c.)X2
% ========> format: column 1 = nodes
% ========> format: column 2 = value

global watchv
xnode=nodesxyz(:,1);
ynode=nodesxyz(:,2);
n1=elemnode(:,1);
n2=elemnode(:,2);
n3=elemnode(:,3);
alphax=elemnode(:,4);
alphay=elemnode(:,5);
beta=elemnode(:,6);
fval=elemnode(:,7);
bnode=boundrnode(:,1);
bval=boundrnode(:,2);

%remove duplicate in boundary nodes
[idxBn,IB,JB]=unique(bnode);
bnode=idxBn;
retBnval=bval(IB);
bval=retBnval;

%initialize lengths
Nelem=length(n1);
Nnode=length(xnode);
Nbound=length(bnode);

%initialize variables
%to store global node at each element
gnode=zeros(3,1);
%to store b^e_i & c^e_i values in computation of Kelem
bei=zeros(3,1);
cei=zeros(3,1);
%to store global K matrix nodexnode
gKmat=sparse([],[],[],Nnode,Nnode);
%to store global B vector
gBvec=zeros(Nnode,1);

for elem_loop=1:Nelem
  %get global node of the current element
  gnode(1)=n1(elem_loop);
  gnode(2)=n2(elem_loop);
  gnode(3)=n3(elem_loop);
  %compute coordinates
  x1=xnode(gnode(1));		    
  x2=xnode(gnode(2));		    
  x3=xnode(gnode(3));		    
  y1=ynode(gnode(1));
  y2=ynode(gnode(2));		    
  y3=ynode(gnode(3));		    

  %compute b^e_i and c^e_i
  bei(1)=y2-y3;
  bei(2)=y3-y1;
  bei(3)=y1-y2;
  cei(1)=x3-x2;
  cei(2)=x1-x3;
  cei(3)=x2-x1;
  %compute area
  area=0.5*abs((bei(1)*cei(2)-bei(2)*cei(1)));
  onePdelta_ij=1+eye(3);
  Ksmall=(alphax(elem_loop).*(bei*bei')+...
	  alphay(elem_loop).*(cei*cei'))./(4.*area)+...
      beta(elem_loop).*onePdelta_ij.*area./12.0;
  gKmat(gnode,gnode)+=Ksmall;
  gBvec(gnode)+=fval(elem_loop).*area./3;
  end
  
end

%set boundary conditions
idx=(1:1:Nnode)';
idxComp=setxor(idx,idxBn);
K=gKmat(idxComp,idxComp);
Kbc=gKmat(idxComp,idxBn);
B=gBvec(idxComp)-Kbc*bval;

%solve system of equation
y=K\B;
